/*
 [auto_generated]
 boost/numeric/odeint/algebra/default_operations.hpp

 [begin_description]
 Default operations. They work with the default numerical types, like float, double, complex< double> ...
 [end_description]

 Copyright 2009-2011 Karsten Ahnert
 Copyright 2009-2011 Mario Mulansky

 Distributed under the Boost Software License, Version 1.0.
 (See accompanying file LICENSE_1_0.txt or
 copy at http://www.boost.org/LICENSE_1_0.txt)
 */

#ifndef BOOST_NUMERIC_ODEINT_ALGEBRA_DEFAULT_OPERATIONS_HPP_INCLUDED
#define BOOST_NUMERIC_ODEINT_ALGEBRA_DEFAULT_OPERATIONS_HPP_INCLUDED

#include <algorithm>

#include <boost/config.hpp>
#include <boost/array.hpp>

#include <boost/numeric/odeint/util/unit_helper.hpp>

namespace boost {
namespace numeric {
namespace odeint {

/*
 * Notes:
 *
 * * the results structs are needed in order to work with fusion_algebra
 */
struct default_operations {

  template <class Fac1 = double>
  struct scale {
    const Fac1 m_alpha1;

    scale(Fac1 alpha1) : m_alpha1(alpha1) {
    }

    template <class T1>
    void operator()(T1& t1) const {
      t1 *= m_alpha1;
    }

    typedef void result_type;
  };

  template <class Fac1 = double>
  struct scale_sum1 {
    const Fac1 m_alpha1;

    scale_sum1(Fac1 alpha1) : m_alpha1(alpha1) {
    }

    template <class T1, class T2>
    void operator()(T1& t1, const T2& t2) const {
      t1 = m_alpha1 * t2;
    }

    typedef void result_type;
  };

  template <class Fac1 = double, class Fac2 = Fac1>
  struct scale_sum2 {
    const Fac1 m_alpha1;
    const Fac2 m_alpha2;

    scale_sum2(Fac1 alpha1, Fac2 alpha2) : m_alpha1(alpha1), m_alpha2(alpha2) {
    }

    template <class T1, class T2, class T3>
    void operator()(T1& t1, const T2& t2, const T3& t3) const {
      t1 = m_alpha1 * t2 + m_alpha2 * t3;
    }

    typedef void result_type;
  };

  template <class Fac1 = double, class Fac2 = Fac1, class Fac3 = Fac2>
  struct scale_sum3 {
    const Fac1 m_alpha1;
    const Fac2 m_alpha2;
    const Fac3 m_alpha3;

    scale_sum3(Fac1 alpha1, Fac2 alpha2, Fac3 alpha3)
      : m_alpha1(alpha1), m_alpha2(alpha2), m_alpha3(alpha3) {
    }

    template <class T1, class T2, class T3, class T4>
    void operator()(T1& t1, const T2& t2, const T3& t3, const T4& t4) const {
      t1 = m_alpha1 * t2 + m_alpha2 * t3 + m_alpha3 * t4;
    }

    typedef void result_type;
  };

  template <class Fac1 = double, class Fac2 = Fac1, class Fac3 = Fac2, class Fac4 = Fac3>
  struct scale_sum4 {
    const Fac1 m_alpha1;
    const Fac2 m_alpha2;
    const Fac3 m_alpha3;
    const Fac4 m_alpha4;

    scale_sum4(Fac1 alpha1, Fac2 alpha2, Fac3 alpha3, Fac4 alpha4)
      : m_alpha1(alpha1), m_alpha2(alpha2), m_alpha3(alpha3), m_alpha4(alpha4) {
    }

    template <class T1, class T2, class T3, class T4, class T5>
    void operator()(T1& t1, const T2& t2, const T3& t3, const T4& t4, const T5& t5) const {
      t1 = m_alpha1 * t2 + m_alpha2 * t3 + m_alpha3 * t4 + m_alpha4 * t5;
    }

    typedef void result_type;
  };

  template <class Fac1 = double, class Fac2 = Fac1, class Fac3 = Fac2, class Fac4 = Fac3,
            class Fac5 = Fac4>
  struct scale_sum5 {
    const Fac1 m_alpha1;
    const Fac2 m_alpha2;
    const Fac3 m_alpha3;
    const Fac4 m_alpha4;
    const Fac5 m_alpha5;

    scale_sum5(Fac1 alpha1, Fac2 alpha2, Fac3 alpha3, Fac4 alpha4, Fac5 alpha5)
      : m_alpha1(alpha1), m_alpha2(alpha2), m_alpha3(alpha3), m_alpha4(alpha4), m_alpha5(alpha5) {
    }

    template <class T1, class T2, class T3, class T4, class T5, class T6>
    void operator()(T1& t1, const T2& t2, const T3& t3, const T4& t4, const T5& t5, const T6& t6) const {
      t1 = m_alpha1 * t2 + m_alpha2 * t3 + m_alpha3 * t4 + m_alpha4 * t5 + m_alpha5 * t6;
    }

    typedef void result_type;
  };

  template <class Fac1 = double, class Fac2 = Fac1, class Fac3 = Fac2, class Fac4 = Fac3,
            class Fac5 = Fac4, class Fac6 = Fac5>
  struct scale_sum6 {
    const Fac1 m_alpha1;
    const Fac2 m_alpha2;
    const Fac3 m_alpha3;
    const Fac4 m_alpha4;
    const Fac5 m_alpha5;
    const Fac6 m_alpha6;

    scale_sum6(Fac1 alpha1, Fac2 alpha2, Fac3 alpha3, Fac4 alpha4, Fac5 alpha5, Fac6 alpha6)
      : m_alpha1(alpha1)
      , m_alpha2(alpha2)
      , m_alpha3(alpha3)
      , m_alpha4(alpha4)
      , m_alpha5(alpha5)
      , m_alpha6(alpha6) {
    }

    template <class T1, class T2, class T3, class T4, class T5, class T6, class T7>
    void operator()(T1& t1, const T2& t2, const T3& t3, const T4& t4, const T5& t5, const T6& t6,
                    const T7& t7) const {
      t1 = m_alpha1 * t2 + m_alpha2 * t3 + m_alpha3 * t4 + m_alpha4 * t5 + m_alpha5 * t6 + m_alpha6 * t7;
    }

    typedef void result_type;
  };

  template <class Fac1 = double, class Fac2 = Fac1, class Fac3 = Fac2, class Fac4 = Fac3,
            class Fac5 = Fac4, class Fac6 = Fac5, class Fac7 = Fac6>
  struct scale_sum7 {
    const Fac1 m_alpha1;
    const Fac2 m_alpha2;
    const Fac3 m_alpha3;
    const Fac4 m_alpha4;
    const Fac5 m_alpha5;
    const Fac6 m_alpha6;
    const Fac7 m_alpha7;

    scale_sum7(Fac1 alpha1, Fac2 alpha2, Fac3 alpha3, Fac4 alpha4, Fac5 alpha5, Fac6 alpha6, Fac7 alpha7)
      : m_alpha1(alpha1)
      , m_alpha2(alpha2)
      , m_alpha3(alpha3)
      , m_alpha4(alpha4)
      , m_alpha5(alpha5)
      , m_alpha6(alpha6)
      , m_alpha7(alpha7) {
    }

    template <class T1, class T2, class T3, class T4, class T5, class T6, class T7, class T8>
    void operator()(T1& t1, const T2& t2, const T3& t3, const T4& t4, const T5& t5, const T6& t6,
                    const T7& t7, const T8& t8) const {
      t1 = m_alpha1 * t2 + m_alpha2 * t3 + m_alpha3 * t4 + m_alpha4 * t5 + m_alpha5 * t6 +
          m_alpha6 * t7 + m_alpha7 * t8;
    }

    typedef void result_type;
  };

  template <class Fac1 = double, class Fac2 = Fac1, class Fac3 = Fac2, class Fac4 = Fac3,
            class Fac5 = Fac4, class Fac6 = Fac5, class Fac7 = Fac6, class Fac8 = Fac7>
  struct scale_sum8 {
    const Fac1 m_alpha1;
    const Fac2 m_alpha2;
    const Fac3 m_alpha3;
    const Fac4 m_alpha4;
    const Fac5 m_alpha5;
    const Fac6 m_alpha6;
    const Fac7 m_alpha7;
    const Fac8 m_alpha8;

    scale_sum8(Fac1 alpha1, Fac2 alpha2, Fac3 alpha3, Fac4 alpha4, Fac5 alpha5, Fac6 alpha6, Fac7 alpha7,
               Fac8 alpha8)
      : m_alpha1(alpha1)
      , m_alpha2(alpha2)
      , m_alpha3(alpha3)
      , m_alpha4(alpha4)
      , m_alpha5(alpha5)
      , m_alpha6(alpha6)
      , m_alpha7(alpha7)
      , m_alpha8(alpha8) {
    }

    template <class T1, class T2, class T3, class T4, class T5, class T6, class T7, class T8, class T9>
    void operator()(T1& t1, const T2& t2, const T3& t3, const T4& t4, const T5& t5, const T6& t6,
                    const T7& t7, const T8& t8, const T9& t9) const {
      t1 = m_alpha1 * t2 + m_alpha2 * t3 + m_alpha3 * t4 + m_alpha4 * t5 + m_alpha5 * t6 +
          m_alpha6 * t7 + m_alpha7 * t8 + m_alpha8 * t9;
    }

    typedef void result_type;
  };

  template <class Fac1 = double, class Fac2 = Fac1, class Fac3 = Fac2, class Fac4 = Fac3,
            class Fac5 = Fac4, class Fac6 = Fac5, class Fac7 = Fac6, class Fac8 = Fac7,
            class Fac9 = Fac8>
  struct scale_sum9 {
    const Fac1 m_alpha1;
    const Fac2 m_alpha2;
    const Fac3 m_alpha3;
    const Fac4 m_alpha4;
    const Fac5 m_alpha5;
    const Fac6 m_alpha6;
    const Fac7 m_alpha7;
    const Fac8 m_alpha8;
    const Fac9 m_alpha9;

    scale_sum9(Fac1 alpha1, Fac2 alpha2, Fac3 alpha3, Fac4 alpha4, Fac5 alpha5, Fac6 alpha6, Fac7 alpha7,
               Fac8 alpha8, Fac9 alpha9)
      : m_alpha1(alpha1)
      , m_alpha2(alpha2)
      , m_alpha3(alpha3)
      , m_alpha4(alpha4)
      , m_alpha5(alpha5)
      , m_alpha6(alpha6)
      , m_alpha7(alpha7)
      , m_alpha8(alpha8)
      , m_alpha9(alpha9) {
    }

    template <class T1, class T2, class T3, class T4, class T5, class T6, class T7, class T8, class T9,
              class T10>
    void operator()(T1& t1, const T2& t2, const T3& t3, const T4& t4, const T5& t5, const T6& t6,
                    const T7& t7, const T8& t8, const T9& t9, const T10& t10) const {
      t1 = m_alpha1 * t2 + m_alpha2 * t3 + m_alpha3 * t4 + m_alpha4 * t5 + m_alpha5 * t6 +
          m_alpha6 * t7 + m_alpha7 * t8 + m_alpha8 * t9 + m_alpha9 * t10;
    }

    typedef void result_type;
  };

  template <class Fac1 = double, class Fac2 = Fac1, class Fac3 = Fac2, class Fac4 = Fac3,
            class Fac5 = Fac4, class Fac6 = Fac5, class Fac7 = Fac6, class Fac8 = Fac7,
            class Fac9 = Fac8, class Fac10 = Fac9>
  struct scale_sum10 {
    const Fac1 m_alpha1;
    const Fac2 m_alpha2;
    const Fac3 m_alpha3;
    const Fac4 m_alpha4;
    const Fac5 m_alpha5;
    const Fac6 m_alpha6;
    const Fac7 m_alpha7;
    const Fac8 m_alpha8;
    const Fac9 m_alpha9;
    const Fac10 m_alpha10;

    scale_sum10(Fac1 alpha1, Fac2 alpha2, Fac3 alpha3, Fac4 alpha4, Fac5 alpha5, Fac6 alpha6,
                Fac7 alpha7, Fac8 alpha8, Fac9 alpha9, Fac10 alpha10)
      : m_alpha1(alpha1)
      , m_alpha2(alpha2)
      , m_alpha3(alpha3)
      , m_alpha4(alpha4)
      , m_alpha5(alpha5)
      , m_alpha6(alpha6)
      , m_alpha7(alpha7)
      , m_alpha8(alpha8)
      , m_alpha9(alpha9)
      , m_alpha10(alpha10) {
    }

    template <class T1, class T2, class T3, class T4, class T5, class T6, class T7, class T8, class T9,
              class T10, class T11>
    void operator()(T1& t1, const T2& t2, const T3& t3, const T4& t4, const T5& t5, const T6& t6,
                    const T7& t7, const T8& t8, const T9& t9, const T10& t10, const T11& t11) const {
      t1 = m_alpha1 * t2 + m_alpha2 * t3 + m_alpha3 * t4 + m_alpha4 * t5 + m_alpha5 * t6 +
          m_alpha6 * t7 + m_alpha7 * t8 + m_alpha8 * t9 + m_alpha9 * t10 + m_alpha10 * t11;
    }

    typedef void result_type;
  };

  template <class Fac1 = double, class Fac2 = Fac1, class Fac3 = Fac2, class Fac4 = Fac3,
            class Fac5 = Fac4, class Fac6 = Fac5, class Fac7 = Fac6, class Fac8 = Fac7,
            class Fac9 = Fac8, class Fac10 = Fac9, class Fac11 = Fac10>
  struct scale_sum11 {
    const Fac1 m_alpha1;
    const Fac2 m_alpha2;
    const Fac3 m_alpha3;
    const Fac4 m_alpha4;
    const Fac5 m_alpha5;
    const Fac6 m_alpha6;
    const Fac7 m_alpha7;
    const Fac8 m_alpha8;
    const Fac9 m_alpha9;
    const Fac10 m_alpha10;
    const Fac11 m_alpha11;

    scale_sum11(Fac1 alpha1, Fac2 alpha2, Fac3 alpha3, Fac4 alpha4, Fac5 alpha5, Fac6 alpha6,
                Fac7 alpha7, Fac8 alpha8, Fac9 alpha9, Fac10 alpha10, Fac11 alpha11)
      : m_alpha1(alpha1)
      , m_alpha2(alpha2)
      , m_alpha3(alpha3)
      , m_alpha4(alpha4)
      , m_alpha5(alpha5)
      , m_alpha6(alpha6)
      , m_alpha7(alpha7)
      , m_alpha8(alpha8)
      , m_alpha9(alpha9)
      , m_alpha10(alpha10)
      , m_alpha11(alpha11) {
    }

    template <class T1, class T2, class T3, class T4, class T5, class T6, class T7, class T8, class T9,
              class T10, class T11, class T12>
    void operator()(T1& t1, const T2& t2, const T3& t3, const T4& t4, const T5& t5, const T6& t6,
                    const T7& t7, const T8& t8, const T9& t9, const T10& t10, const T11& t11,
                    const T12& t12) const {
      t1 = m_alpha1 * t2 + m_alpha2 * t3 + m_alpha3 * t4 + m_alpha4 * t5 + m_alpha5 * t6 +
          m_alpha6 * t7 + m_alpha7 * t8 + m_alpha8 * t9 + m_alpha9 * t10 + m_alpha10 * t11 +
          m_alpha11 * t12;
    }

    typedef void result_type;
  };

  template <class Fac1 = double, class Fac2 = Fac1, class Fac3 = Fac2, class Fac4 = Fac3,
            class Fac5 = Fac4, class Fac6 = Fac5, class Fac7 = Fac6, class Fac8 = Fac7,
            class Fac9 = Fac8, class Fac10 = Fac9, class Fac11 = Fac10, class Fac12 = Fac11>
  struct scale_sum12 {
    const Fac1 m_alpha1;
    const Fac2 m_alpha2;
    const Fac3 m_alpha3;
    const Fac4 m_alpha4;
    const Fac5 m_alpha5;
    const Fac6 m_alpha6;
    const Fac7 m_alpha7;
    const Fac8 m_alpha8;
    const Fac9 m_alpha9;
    const Fac10 m_alpha10;
    const Fac11 m_alpha11;
    const Fac12 m_alpha12;

    scale_sum12(Fac1 alpha1, Fac2 alpha2, Fac3 alpha3, Fac4 alpha4, Fac5 alpha5, Fac6 alpha6,
                Fac7 alpha7, Fac8 alpha8, Fac9 alpha9, Fac10 alpha10, Fac11 alpha11, Fac12 alpha12)
      : m_alpha1(alpha1)
      , m_alpha2(alpha2)
      , m_alpha3(alpha3)
      , m_alpha4(alpha4)
      , m_alpha5(alpha5)
      , m_alpha6(alpha6)
      , m_alpha7(alpha7)
      , m_alpha8(alpha8)
      , m_alpha9(alpha9)
      , m_alpha10(alpha10)
      , m_alpha11(alpha11)
      , m_alpha12(alpha12) {
    }

    template <class T1, class T2, class T3, class T4, class T5, class T6, class T7, class T8, class T9,
              class T10, class T11, class T12, class T13>
    void operator()(T1& t1, const T2& t2, const T3& t3, const T4& t4, const T5& t5, const T6& t6,
                    const T7& t7, const T8& t8, const T9& t9, const T10& t10, const T11& t11,
                    const T12& t12, const T13& t13) const {
      t1 = m_alpha1 * t2 + m_alpha2 * t3 + m_alpha3 * t4 + m_alpha4 * t5 + m_alpha5 * t6 +
          m_alpha6 * t7 + m_alpha7 * t8 + m_alpha8 * t9 + m_alpha9 * t10 + m_alpha10 * t11 +
          m_alpha11 * t12 + m_alpha12 * t13;
    }

    typedef void result_type;
  };

  template <class Fac1 = double, class Fac2 = Fac1, class Fac3 = Fac2, class Fac4 = Fac3,
            class Fac5 = Fac4, class Fac6 = Fac5, class Fac7 = Fac6, class Fac8 = Fac7,
            class Fac9 = Fac8, class Fac10 = Fac9, class Fac11 = Fac10, class Fac12 = Fac11,
            class Fac13 = Fac12>
  struct scale_sum13 {
    const Fac1 m_alpha1;
    const Fac2 m_alpha2;
    const Fac3 m_alpha3;
    const Fac4 m_alpha4;
    const Fac5 m_alpha5;
    const Fac6 m_alpha6;
    const Fac7 m_alpha7;
    const Fac8 m_alpha8;
    const Fac9 m_alpha9;
    const Fac10 m_alpha10;
    const Fac11 m_alpha11;
    const Fac12 m_alpha12;
    const Fac13 m_alpha13;

    scale_sum13(Fac1 alpha1, Fac2 alpha2, Fac3 alpha3, Fac4 alpha4, Fac5 alpha5, Fac6 alpha6,
                Fac7 alpha7, Fac8 alpha8, Fac9 alpha9, Fac10 alpha10, Fac11 alpha11, Fac12 alpha12,
                Fac13 alpha13)
      : m_alpha1(alpha1)
      , m_alpha2(alpha2)
      , m_alpha3(alpha3)
      , m_alpha4(alpha4)
      , m_alpha5(alpha5)
      , m_alpha6(alpha6)
      , m_alpha7(alpha7)
      , m_alpha8(alpha8)
      , m_alpha9(alpha9)
      , m_alpha10(alpha10)
      , m_alpha11(alpha11)
      , m_alpha12(alpha12)
      , m_alpha13(alpha13) {
    }

    template <class T1, class T2, class T3, class T4, class T5, class T6, class T7, class T8, class T9,
              class T10, class T11, class T12, class T13, class T14>
    void operator()(T1& t1, const T2& t2, const T3& t3, const T4& t4, const T5& t5, const T6& t6,
                    const T7& t7, const T8& t8, const T9& t9, const T10& t10, const T11& t11,
                    const T12& t12, const T13& t13, const T14& t14) const {
      t1 = m_alpha1 * t2 + m_alpha2 * t3 + m_alpha3 * t4 + m_alpha4 * t5 + m_alpha5 * t6 +
          m_alpha6 * t7 + m_alpha7 * t8 + m_alpha8 * t9 + m_alpha9 * t10 + m_alpha10 * t11 +
          m_alpha11 * t12 + m_alpha12 * t13 + m_alpha13 * t14;
    }

    typedef void result_type;
  };

  template <class Fac1 = double, class Fac2 = Fac1, class Fac3 = Fac2, class Fac4 = Fac3,
            class Fac5 = Fac4, class Fac6 = Fac5, class Fac7 = Fac6, class Fac8 = Fac7,
            class Fac9 = Fac8, class Fac10 = Fac9, class Fac11 = Fac10, class Fac12 = Fac11,
            class Fac13 = Fac12, class Fac14 = Fac13>
  struct scale_sum14 {
    const Fac1 m_alpha1;
    const Fac2 m_alpha2;
    const Fac3 m_alpha3;
    const Fac4 m_alpha4;
    const Fac5 m_alpha5;
    const Fac6 m_alpha6;
    const Fac7 m_alpha7;
    const Fac8 m_alpha8;
    const Fac9 m_alpha9;
    const Fac10 m_alpha10;
    const Fac11 m_alpha11;
    const Fac12 m_alpha12;
    const Fac13 m_alpha13;
    const Fac14 m_alpha14;

    scale_sum14(Fac1 alpha1, Fac2 alpha2, Fac3 alpha3, Fac4 alpha4, Fac5 alpha5, Fac6 alpha6,
                Fac7 alpha7, Fac8 alpha8, Fac9 alpha9, Fac10 alpha10, Fac11 alpha11, Fac12 alpha12,
                Fac13 alpha13, Fac14 alpha14)
      : m_alpha1(alpha1)
      , m_alpha2(alpha2)
      , m_alpha3(alpha3)
      , m_alpha4(alpha4)
      , m_alpha5(alpha5)
      , m_alpha6(alpha6)
      , m_alpha7(alpha7)
      , m_alpha8(alpha8)
      , m_alpha9(alpha9)
      , m_alpha10(alpha10)
      , m_alpha11(alpha11)
      , m_alpha12(alpha12)
      , m_alpha13(alpha13)
      , m_alpha14(alpha14) {
    }

    template <class T1, class T2, class T3, class T4, class T5, class T6, class T7, class T8, class T9,
              class T10, class T11, class T12, class T13, class T14, class T15>
    void operator()(T1& t1, const T2& t2, const T3& t3, const T4& t4, const T5& t5, const T6& t6,
                    const T7& t7, const T8& t8, const T9& t9, const T10& t10, const T11& t11,
                    const T12& t12, const T13& t13, const T14& t14, const T15& t15) const {
      t1 = m_alpha1 * t2 + m_alpha2 * t3 + m_alpha3 * t4 + m_alpha4 * t5 + m_alpha5 * t6 +
          m_alpha6 * t7 + m_alpha7 * t8 + m_alpha8 * t9 + m_alpha9 * t10 + m_alpha10 * t11 +
          m_alpha11 * t12 + m_alpha12 * t13 + m_alpha13 * t14 + m_alpha14 * t15;
    }

    typedef void result_type;
  };

  template <class Fac1 = double, class Fac2 = Fac1>
  struct scale_sum_swap2 {
    const Fac1 m_alpha1;
    const Fac2 m_alpha2;

    scale_sum_swap2(Fac1 alpha1, Fac2 alpha2) : m_alpha1(alpha1), m_alpha2(alpha2) {
    }

    template <class T1, class T2, class T3>
    void operator()(T1& t1, T2& t2, const T3& t3) const {
      const T1 tmp(t1);
      t1 = m_alpha1 * t2 + m_alpha2 * t3;
      t2 = tmp;
    }

    typedef void result_type;
  };

  /*
   * for usage in for_each2
   *
   * Works with boost::units by eliminating the unit
   */
  template <class Fac1 = double>
  struct rel_error {
    const Fac1 m_eps_abs, m_eps_rel, m_a_x, m_a_dxdt;

    rel_error(Fac1 eps_abs, Fac1 eps_rel, Fac1 a_x, Fac1 a_dxdt)
      : m_eps_abs(eps_abs), m_eps_rel(eps_rel), m_a_x(a_x), m_a_dxdt(a_dxdt) {
    }

    template <class T1, class T2, class T3>
    void operator()(T3& t3, const T1& t1, const T2& t2) const {
      using std::abs;
      set_unit_value(t3, abs(get_unit_value(t3)) / (m_eps_abs +
                                                    m_eps_rel * (m_a_x * abs(get_unit_value(t1)) +
                                                                 m_a_dxdt * abs(get_unit_value(t2)))));
    }

    typedef void result_type;
  };

  /*
   * for usage in for_each3
   *
   * used in the controller for the rosenbrock4 method
   *
   * Works with boost::units by eliminating the unit
   */
  template <class Fac1 = double>
  struct default_rel_error {
    const Fac1 m_eps_abs, m_eps_rel;

    default_rel_error(Fac1 eps_abs, Fac1 eps_rel) : m_eps_abs(eps_abs), m_eps_rel(eps_rel) {
    }

    /*
     * xerr = xerr / ( eps_abs + eps_rel * max( x , x_old ) )
     */
    template <class T1, class T2, class T3>
    void operator()(T3& t3, const T1& t1, const T2& t2) const {
      BOOST_USING_STD_MAX();
      using std::abs;
      Fac1 x1 = abs(get_unit_value(t1)), x2 = abs(get_unit_value(t2));
      set_unit_value(t3, abs(get_unit_value(t3)) /
                         (m_eps_abs + m_eps_rel * max BOOST_PREVENT_MACRO_SUBSTITUTION(x1, x2)));
    }

    typedef void result_type;
  };

  /*
   * for usage in reduce
   */

  template <class Value>
  struct maximum {
    template <class Fac1, class Fac2>
    Value operator()(Fac1 t1, const Fac2 t2) const {
      using std::abs;
      Value a1 = abs(get_unit_value(t1)), a2 = abs(get_unit_value(t2));
      return (a1 < a2) ? a2 : a1;
    }

    typedef Value result_type;
  };

  template <class Fac1 = double>
  struct rel_error_max {
    const Fac1 m_eps_abs, m_eps_rel;

    rel_error_max(Fac1 eps_abs, Fac1 eps_rel) : m_eps_abs(eps_abs), m_eps_rel(eps_rel) {
    }

    template <class Res, class T1, class T2, class T3>
    Res operator()(Res r, const T1& x_old, const T2& x, const T3& x_err) {
      BOOST_USING_STD_MAX();
      using std::abs;
      Res tmp = abs(get_unit_value(x_err)) /
          (m_eps_abs + m_eps_rel * max BOOST_PREVENT_MACRO_SUBSTITUTION(abs(x_old), abs(x)));
      return max BOOST_PREVENT_MACRO_SUBSTITUTION(r, tmp);
    }
  };

  template <class Fac1 = double>
  struct rel_error_max2 {
    const Fac1 m_eps_abs, m_eps_rel, m_a_x, m_a_dxdt;

    rel_error_max2(Fac1 eps_abs, Fac1 eps_rel, Fac1 a_x, Fac1 a_dxdt)
      : m_eps_abs(eps_abs), m_eps_rel(eps_rel), m_a_x(a_x), m_a_dxdt(a_dxdt) {
    }

    template <class Res, class T1, class T2, class T3, class T4>
    Res operator()(Res r, const T1& x_old, const T2& x, const T3& dxdt_old, const T4& x_err) {
      BOOST_USING_STD_MAX();
      using std::abs;
      Res tmp = abs(get_unit_value(x_err)) /
          (m_eps_abs +
           m_eps_rel * (m_a_x * abs(get_unit_value(x_old)) + m_a_dxdt * abs(get_unit_value(dxdt_old))));
      return max BOOST_PREVENT_MACRO_SUBSTITUTION(r, tmp);
    }
  };

  template <class Fac1 = double>
  struct rel_error_l2 {
    const Fac1 m_eps_abs, m_eps_rel;

    rel_error_l2(Fac1 eps_abs, Fac1 eps_rel) : m_eps_abs(eps_abs), m_eps_rel(eps_rel) {
    }

    template <class Res, class T1, class T2, class T3>
    Res operator()(Res r, const T1& x_old, const T2& x, const T3& x_err) {
      BOOST_USING_STD_MAX();
      using std::abs;
      Res tmp = abs(get_unit_value(x_err)) /
          (m_eps_abs + m_eps_rel * max BOOST_PREVENT_MACRO_SUBSTITUTION(abs(x_old), abs(x)));
      return r + tmp * tmp;
    }
  };

  template <class Fac1 = double>
  struct rel_error_l2_2 {
    const Fac1 m_eps_abs, m_eps_rel, m_a_x, m_a_dxdt;

    rel_error_l2_2(Fac1 eps_abs, Fac1 eps_rel, Fac1 a_x, Fac1 a_dxdt)
      : m_eps_abs(eps_abs), m_eps_rel(eps_rel), m_a_x(a_x), m_a_dxdt(a_dxdt) {
    }

    template <class Res, class T1, class T2, class T3, class T4>
    Res operator()(Res r, const T1& x_old, const T2& x, const T3& dxdt_old, const T4& x_err) {
      using std::abs;
      Res tmp = abs(get_unit_value(x_err)) /
          (m_eps_abs +
           m_eps_rel * (m_a_x * abs(get_unit_value(x_old)) + m_a_dxdt * abs(get_unit_value(dxdt_old))));
      return r + tmp * tmp;
    }
  };
};

}  // odeint
}  // numeric
}  // boost

#endif  // BOOST_NUMERIC_ODEINT_ALGEBRA_DEFAULT_OPERATIONS_HPP_INCLUDED
